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indexes (valid from January 1997) 
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Abstract. Phase mixing of a three dimensional magnetohydrodynamic (MHD) pulse is studied in the compressive, 
three-dimensional (without an ignorable coordinate) regime. It is shown that the efficiency of decay of an Alfvenic 
part of a compressible MHD pulse is related linearly to the degree of localization of the pulse in the homogeneous 
transverse direction. In the developed stage of phase mixing (for large times), coupling to its compressive part 
does not alter the power-law decay of an Alfvenic part of a compressible MHD pulse. The same applies to the 
dependence upon the resistivity of the Alfvenic part of the pulse. All this implies that the dynamics of Alfven 
waves can still be qualitatively understood in terms of the previous 2.5D models. Thus, the phase mixing remains 
a relevant paradigm for the coronal heating applications in the realistic 3D geometry and compressive plasma. 

Key words. Magnetohydrodynamics(MHD)- waves - Sun: activity - Sun: Corona 



1. Introduction 

The phase mixing of incompressible Alfven waves 
ijHevvaerts fc Priest I Il983ii has been intensively studied 
as a possible mechanism of coronal heating in open mag- 
netic structures of the polar regions of the Sun. Also, there 
are recent observational indicatio ns that the same mech- 
anism operates in coronal loops lOfman fc Aschwanden I 
ll2002t) . In a compressible plasma, phase mixing of lin- 
early polarized plane Alfven waves may provide enhanced 
heating through the nonlinear generation of fast magne- 
toacoustic waves. This p henomen on has been in tensively 
studi e d bvlMalara et al. I (d996'l. Nakariakov et al. 1199'^ . 



Il998|) . iBotha et al. I l|200Cl|) . iTsiklauri et al. I l|200lt 12^' 

both analytically and numerically. 

There are three main motivations for further devel- 
op ment of this work: T he first, as was recently shown 
by iHood et al. I l)2002a|) . is that phase mixing of local- 
ized single Alfven pulses results in a slower, power-law 
damping as opposed to the standard oc exp(— t'^) on e 
for harmonic Alfven waves l)Hevvaerts fc Priest I Il98^ . 
This suggests that localized Alfvenic perturbations can 
transport energy higher into th e corona than harmonic 
ones. Also, "Hoo d et al. I l)2002bfl showed that this power- 
law damping is replaced by even faster power-law decay 
as the number of pulses is increased. Such perturbations 
could be generated e.g. by transie nt events such as sola r 
flares, coronal mass ejections, etc (iRoussev et al. lE00l|) . 
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The second follows from the recent demonstration by 
Irsiklauri fc Nakariakov ( 2002) that in the linear regime 
and in 3D geometry, generation of compressive perturba- 
tions significantly modifies the dynamics of the Alfvenic 
part of a MHD pulse. In fact, while solving the propaga- 
tion part of the problem, i.e. considering an ideal plasma 
limit, it has been shown that an initially localized 3D 
Alfven pulse interacting with a ID transverse inhomogene- 
ity evolves to a MHD pulse with a significant compressible 
component. This 3D MHD pulse is subject to phase mix- 
ing. Therefore, as a natural step forward, in this study we 
include a finite plasma resistivity in order to investigate 
quantitatively dissipation of a 3D MHD pulse via phase 
mixing. The third follows from Parker's (1991) earlier sug- 
gestion that in order for the phase mixing mechanism to 
work, it requires the presence of a third ignorable coor- 
dinate. Although it has been justly observed that having 
an ignorable coordinate seems unlikely in a filamentary 
corona, and, also, in the absence of an ignorable coordi- 
nate, the "Alfven-type waves" are coupled to the "fast 
magnetosonic-type modes", we believe it was not justi- 
fied to dismiss the importance of phase mixing without a 
rigorous, fully three dimensional study. 

2. The model 

In our model we use the MHD equations 

p dtV + p{V ■ V)V = -Vp - {At:)-^B x curlS, (1) 

dtB = curl(F X B) -t- ??AB, (2) 
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Fig. 1. Time evolution of max{By{x = 0,y,z,t) (solid 
lines) and max{Vy{x = 0,y,z,t) (dotted lines). Note the 
plot is log- normal. Here, ay = 1.0. The thick lines rep- 
resent simulation results with the doubled resolution in 
y-direction. Note progressive convergence of the solutions 
with the increase of simulation time. 

dtP + V ■Vp + 7pV • V = (7 - l)ri{cuT\Bf, (3) 

where B is the magnetic field, V is the plasma velocity, p 
is the plasma mass density, and p is the plasma thermal 
pressure. In what follows we use 5/3 for the value of 7. 
Also, here rj denotes plasma resistivity. 

We solve equations Q-© in Cartesian coordinates 
{x,y,z). Note that we solve a fully 3D problem retaining 
variation in the y-direction, i.e. {d/dy ^ 0). A uniform 
magnetic field Bq is in the z-direction. The plasma config- 
uration has a one-dimensional inhomogeneity in the equi- 
librium density poix) and temperature Tq{x), while the 
unperturbed thermal pressure, po, is taken to be constant 
everywhere. 

3. Numerical Results 

Evolution of a purely Alfvenic perturbation allows a full 
analytical treatment and this is presented in Appendix 
A. However, if an initial Alfvenic perturbation depends 
on the y-coordinatc, contrary to the case studied in 
the Appendix A, the linearized version of the system of 
equations Q-® no longer allows decoupling of Alfvenic 
and compressive perturbations. In fact, all three modes, 
Alfven, fast and slow magnetosonic, are inter-coupled. 
Thus, they cannot be separated. Because of the complexity 
of the problem in this case, we resort to direct numerical 
simulations of Eqs.l Ul-lEt using the numerical code LareSd 
(jArber et al. Il200l|) . We employ the following initial con- 
ditions: 




\z\<0.1, (4) 
elsewhere 



Fig. 2. The developed stage of phase mixing: time evo- 
lution of max{By{x = 0,y,z,t) for large times (log-log 
plot). The thick solid line represents the analytic solution 
fEa. (|A.10|) with an additional pre-factor 0.5, to account 
for the split of the initial pulse into two D'Alambert solu- 
tions), while the thin line is our numerical simulation re- 
sult (for ay — 0). Note that the deviation between the two 
progressively vanishes. The dotted, dashed, dash-dotted 
and dash-triple-dotted curves correspond to the numeri- 
cal solutions for ay = 0.5, 1.0, 1.5, 2.0 respectively. 



The rest of the perturbations at t = are set to zero 
(including By, in this way we guarantee fulfillment of 
divB = 0, at t = 0). In Eq.Q A is the amphtude of 
the wave, Ca{x) is the local Alfven speed (see Appendix 
A) and ay is a free parameter that controls the strength 
of the gradient in the y direction of the initial pertur- 
bation. In effect, ay controls the strength of the coupling 
between the Alfvenic and compressive (fast and slow mag- 
netosonic) components of the MHD pulse. This particu- 
lar choice of y-dependence in Eq.Q is dictated by a de- 
sire to have the maximum of the absolute value of dyVy, 
which is responsible for the coupling to the compressive 
modes, at the middle of the simulation cube, where the 
background plasma inhomogeneity is the strongest. Since 
LareSd is a fully non-linear numerical code, we simply set 
small (linear) initial amplitude of the perturbation (given 
by Eq.0)) as A = 0.0001, while the output results are 
normalized to A. The simulation cube size was set by the 
limits —11.0 < x,y,z < 11.0. Boundary conditions used 
in all our simulations are zero-gradient in all three spa- 
tial dimensions. In fact, our strategy was to forbid any 
generated wave fronts to reach the boundaries to avoid 
contamination of the solution. This of course severely re- 
stricts simulation end time. However, with the computa- 
tional resources at our disposal, it was still possible to 
reach the developed stage of phase mixing. In all our nu- 
merical simulations we fix the plasma /3 at 10~^ (which is 
a plausible value for solar corona) and dimensionless re- 
sistivity (see Appen dix A), ii, at 5 x 10~ ^ (which is the 
same value used by iHood et al. I l|2002al) ) unless stated 
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Fig. 3. EfRciency of decay of the Alfvenic part of MHD 

pulse, max{By{x — 0,y, z)/max{By{x = 0,y,z,ay — 
0), as a function of coupling to the compressive modes 
strength parameter, aj, at t = 10. 



CQ 




Fig. 4. Dependence of max{By{x = 0,y, z,t = 10) upon 
resistivity rj. The solid line represents the analytic solution 
(Eg. ljA.lOp with an additional pre-factor 0.5, to account 
for the split of the initial pulse into two D'Alambert so- 
lutions and a further pre-factor 0.7188, which comes from 
the middle point of Fig. (|3Jl), while open diamonds are 
numerical simulation results. Here, a,, = 1.0. 



otherwise. In Fig. iQ} we track the evolution of the max- 
imal value (over the whole simulation cube at a given 
time) of By and Vy at the middle of the density inho- 
mogeneity at a: = 0, where phase mixing significantly en- 
hances the dissipation. Note that max{\Vy{x — 0,y,z,t)\) 
at t = is not unity as it is normalized to A, not to 
max{\Vy{x = 0,y,z,t = 0)|). We gather from this graph 
that although By was initially zero, it is rapidly gener- 
ated (on the time scale of ~ 0.05, which is the ratio of 
length-scale of the pulse I/IOtt to Ca{0) — 0.6), and sub- 
sequently it decays due to phase mixing. Most of our nu- 
merical runs were done with 200 x 72 x 400 resolution 
(with a non- uniform grid), and we checked convergence of 



our results by doubling the number of grid points in the 
y-direction. Coincidence of the two runs was satisfactory 
and, in fact, it progressively improves with the increase of 
the simulation time. Although we were unable to double 
the resolution in x and z, proof of the convergence of our 
results comes from the fact (see Fig. (|2l) that the analytic 
result (given by Ea. ljA.10|l with an additional pre-factor 
0.5 to account for the splitting of initial V„ perturbation 
into two D'Alambert solutions I Tsiklauri fc Nakariakov I 



I2OO2I) ) and the numerical solution in the case ay = 
(when there is a complete decoupling of the Alfven and 
the compressive modes) coincide perfectly for large ts. In 
fact, at t = 9.998143, the absolute deviation between the 
two is (0.077836 - 0.077785)/0.077785 = 0.0006, which 
is a remarkably good match for a fully three dimensional 
simulation. Figure (O demonstrates that: (A) an increase 
in the coupling strength parameter, ay, results in an over- 
all decrease of the Alfvenic part (By) of the MHD pulse, 
as expected, and (B) At large times, all lines on the log- 
log plot tend to become parallel. Therefore, the ampli- 
tude of the Alfvenic component of the MHD pulse de- 
cays with the same power-law {By oc t^^^^) as in the 
case of a localized, pure, Alfven wave (which occurs only 
when ay = 0), i.e. coupling to the compressive waves does 
not affect the decay power-law of the Alfvenic part of the 
MHD pulse. In order to investigate and quantify point 
(A) further, in Fig. Q we plot the eSiciency of decay of 
the Alfvenic part of the MHD pulse, which we define as 
max(By(x = 0,y, z)/max{By{x = 0,y,z,ay = 0), as a 
function of coupling to the compressive modes strength 
parameter, ay at t = 10. Open diamonds represent nu- 
merical simulation results, while a straight line is the best 
fit, 1.0 — 0.25403cky. We conclude that the efficiency of 
decay depends linearly upon the degree of localization of 
the pulse in the homogeneous transverse direction, i.e. the 
coupling strength parameter, ay. As far as point (B) is 
concerned, apart from the visual manifestation that at 
large times all lines tend to become parallel to each other, 
an actual fitting procedure also confirms that, indeed, am- 
plitude of the Alfvenic part of the MHD pulse decays as 
a power-law (cx t~^/^) within a few percent error to the 
power-law index. This is not a trivial result. From this 
we reach the conclusion that the dynamics (phase mixing) 
of weakly non-plane Alfven waves can .still be qualitatively 
understood in terms of the previous 2. 5D models. A certain 
level of understanding of this profound conclusion can be 
achieved by analyzing a linearized version of the system of 
equations Q-© for /3 = 0. In this case, coupling between 
the Alfven {By) and fast magnetosonic mode {B^) (the 
slow mode is absent in this approximation) is described 
by the following system of equations: 



[5; 



Cl{x){dl,+dl)] B, 



vdtdl,B^+ C\{x)dlyBy 



[5- 



C\{x){d'^ 



yy 



-dl)]By-- 
vdtdl^By- 



Cl{x)dlyB^ 



(5) 



(6) 
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Fig. 5. Evolution of \Vy{x = Q,y ~ 0, z)\ z-profile in time. 
Thin line (spike in the middle) is at t = 0, while thick and 
the thickest lines depict \Vy{x = 0,y = 0, z)| at t = 5.0 
and 10.0. Here, = 1.0. 



In order to fully understand the dynamics of the 
Alfvenic part of the MHD pulse under our initial con- 
ditions as well as to appreciate fully the extent of the effi- 
ciency of the phase mixing for the coronal plasma heating, 
in Fig. lO we show the evolution of the absolute value of 
the ^-profile in the middle of the plasma inhomogeneity 
of Vy with time. We gather from this graph that the ini- 
tial, very localized, profile splits into two D'Alambert so- 
luti ons propagating in opposite dir ections along the field 
fcf. lTsiklauri fc Nakariakov I l|2002t) V which then dissipate 
rapidly towards a diffuse, Gaussian shape. 

Finally, we would like to reiterate the importance 
of plasma compressibility on the MHD wave dynamics 
in the context of wave-based t heories of coronal heat- 
ing a nd MHD turbulence (cf. IXsiklauri fc Nakariakov I 
l(2flfl2l) '). Our present numerical simulations confirm that 
when a J, 7^ 0, coupling to the Alfvenic part of the MHD 
pulse quite efficiently generates a density perturbation 
(which was initially zero) on a sizeable fraction of the ini- 
tial Alfven amplitude. 



From Eqs.jSj) and lO it follows that the terms responsi- 
ble for the coupling between the Alfven and fast magne- 
tosonic mode are C'^B^.y/iL^Ly) (where L^ y are spatial 
scales in x and y), while the terms responsible for the dis- 
sipation are rjBx^y / {IaLI.) (where tA is a typical Alfven 
time scale). At large is, because of the action of phase 
mixing in the region of density inhomogeneity, dissipa- 
tive terms (oc L~'^) are much greater than the coupling 
terms (oc L~^). Therefore, for i ^ 1 the dynamics of By 
is governed by the dissipation rather than its coupling to 
the compressive fast mode. This, serves as an explanation 
why, at large times, all lines in Fig. (O tend to become 
parallel, i.e. the amplitude of Alfvenic component of the 
MHD pulse decays as the same power-law {By cx t'^/"^) 
as in the case of localized, pure, Alfven wave. 

Also, we investigated the dependence of the Alfvenic 
part of the MHD pulse, which we quantify by max{By{x — 
0,y,z,t = 10), upon resistivity, 77. The results are pre- 
sented in Fig. 01 On the one hand, an accurate analytic 
solution IIA.10|) for By, which scales as oc l/y^, was ob- 
tained for the case of ay — 0. On the other hand, point 
(B), mentioned above, indicates that when ay =/= 0, the 
solution IjA.lOl) should be altered by a further pre-factor 
that can be determined e.g. from Fig. © (in this case, as 
Qfj, = 1, it is 0.7188). Therefore, we can see that this way 
constructed analytic solution gives a tolerably good fit in 
Fig. The deviation in the small 77 end can be easily 
explained by the fact that the solution (jA.lOp was derived 
in the developed phase mixing approximation, which for 
such small values of rj would require numerical runs with 
end-simulation times much larger than t = 10 (the lat- 
ter is not possible with current computational facilities, 
as it would require further increase of the simulation box 
size). Thus, we conclude that the coupling to its compres- 
sive part does not alter the 1 / ^ behavior of the Alfvenic 
part of a compressible MHD pulse. 



4. Conclusions 

In summary, we have established that: 

— Efficiency of the decay of the Alfvenic part of a com- 
pressible MHD pulse is related linearly to the degree of 
localization of the pulse in the homogeneous transverse 
direction, i.e. the coupling strength parameter ay. 

— In the developed stage of phase mixing (for large 
times) , coupling to its compressive part does not alter 
the l/^ behavior of the Alfvenic part of a compress- 
ible MHD pulse. 

— The dynamics of Alfven waves can still be qualitatively 
understood in terms of the previous 2.5D models, e.g., 
iHood et al.1 l|20024 's law {By oc t'^/'^) for spatially 
localized, single pulses, survives in the case of 3D com- 
pressive MHD waves (in the absence of an ignorable 
coordinate). 

— All of the above implies that phase mixing remains a 
relevant paradigm for the coronal heating applications 
in realistic 3D geometry and compressive plasma. 

It is useful to give estimates of typical scales of 
the problem. We choose a typical (normalization) length 
scale of the transverse inhomogeneity L=10 Mm and 
Alfven speed Ca = 1000 km s^^, which are appropri- 
ate for example for coronal plumes. The coronal value 
for shear (Braginskii) viscosity, by which Alfven waves 
damp, is es t imated to be about rj = 1 m^ 
iHood et al. I l)2002aj) . which translates into 10 
dimensionless units). However, as some recent studies 
sugge st (Nakariakov et al. 19991 lOfman fc Aschwanden I 
I2OO2I) the actual value of the shear viscosity is possibly 
much lar ger due to micro- turbule nce, e.g. IQ^-^^^-^ 
s-^ (O fman fc Aschwanden 1 12002|) . Effectively, this im- 
plies that the Alfven waves damp by bulk Braginskii vis- 
cosity rather than the shear one. The latter estimate is 
based on the observations of standing (kink) Alfven waves 
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in the coronal loops. Thus, its use in our estimates is justi- 
fied provided the propagating Alfven waves, studied here, 
are affected by the same dissipation mechanism as the 
standing ones. Using Ea. ljA.QI) (for a; = 0) we can esti- 
mate the characteristic time of damping of the amplitude 
of an Alfven wave (when it reduces, say, e-times) due to 
phase mixing. Thus, we obtain approximately 2745.22 s 
for 7/ 1 m-^ s"^ and 0.16 — 34.56 s for r/ — ^O^-^^-^-^ 
s^^. This translates into 1647.13 Mm (the Alfven speed in 
the middle of inhomogeneity is 600 km s~^) for the damp- 
ing length, for 7? = 1 m^ s"\ and 0.10 - 20.74 Mm for 
r\ = lo9-2±3-5 ni^ s^^. A similar analysis for the harmonic 
waves (now using Ea. ljA.7() and substituting fc = in /a) 
yields a damping length of 609.47 Mm for 77 = 1 m^ s^^ 
and 0.04 - 7.67 Mm for rj = 109-2±3.5 ^2 g-i. These esti- 
mates would be affected if we choose different (probably 
higher) Alfven speeds or if the pulse is wider. Also, note 
that in our numerical runs, the width of the Gaussian 
Alfven pulse in the longitudinal (z) direction was fixed at 
(7 « 0.03, i.e. 0.3 Mm, while in the transverse (y) direction 
the width, l/y/2ay, was varied from -f 00 down to 0.35, i.e. 
from -\-oo down to 3.5 Mm. 
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Appendix A: 

Here, we present an alternative, more rigorou s, way to de- 
rive th e power-law decay of By , suggested by ijHood et al. I 
l2002a|) . If the initial Alfvenic perturbation (perturbing Vy 
and By) does not depend on the y-coordinate, the lin- 
earized version of the system of equations CJ-® allows 
for a complete decoupling of Alfvenic and compressive per- 
turbations. Therefore, in order to investigate phase mixing 
quantitatively, we solve the following initial value problem 
with zero-gradient boundary conditions: 



[5i 



B 



y\t=o 



B,, 



vdtdl^By 



(A.l) 
(A.2) 



Here, Ca{x) is a dimensionless Alfven speed. In our nu- 
merical simulations we have used the following profile 
for Ca{x): Ca{x) = 1 for x < -0.5, Ca{x) = 0.6 + 
0.4cos[7r(a; 0.5)] for x S [-0.5,0.5], and Ca{x) = 0.2 
fo r a: > 0.5 (we employ here the same normalization as 
in lTsiklauri fc Nakariakov] l)2002(l and, -q is normalized to 
the product of ref erence length and reference Alfven speed 
as in lHood et al. 1 (2002a) ). Note, that as in the developed 
stage of phase mixing (for large times) c?^^. ^ dyy,d1^, in 
the forthcoming analytic treatment, we only retain largest 
(over a;-coordinate) derivatives in the Laplacian (while for 
generality we still keep the full Laplacian in our numeri- 
cal simulations, see below). Next, we introduce the follow- 
ing Lagrangian coordinates and slow (dissipation) time 



scale: x = x, — z — CA{x)t, r = et, with e ^ 1. In 
these variables, the leading term on the left hand side 
of Ea. (|A.l|l is —2eCAix)d'^^By, while on the right hand 
side the leading term is —r]CA{x)d^[C'^{x)'^{T/e)^d'^^By]. 
Here, prime denotes a derivative over x. After integration 
over ^ and introduction of yet another auxiliary variable, 
s = r]C'^{x)^Ty{6e^) = rjC'^{x)H^ /6, we obtain the fol- 
lowing (diffusion) equation for By-. 



dgBy = d^^By, 

which can be easily integrated: 



ITS 



r+00 




/ exp 

J —00 


4s 



(A.3) 



By{e,t = 0)dC'.{AA) 



From Ea. (|A.4p it immediately follows that in the devel- 
oped stage of phase mixing, i.e. when t, s — > +00, the 
solution reduces to 



1 /■+°° 
Therefore, provided the integral 

By{e,t^O)d^', 



(A.5) 



(A.6) 



is finite (which is not always so e.g. for the harmonic initial 
conditions, see below). By scales as l/-\/s or equivalently 

Sj, OC77-l/2t-3/2. 

If we substitute a harmonic wave with initial condi- 
tions By{^',t = 0) = exp(ifc^') into Ea. HA.4p . simple inte- 
gration yields 



B„ 



^-■nC'A(xfee /6 ^-ik(z-CA{x)t) 



(A.7) 



thus, we recover the well-known Heyvaerts & Priest's so- 
lution. For the spatially localized (in z) initial conditions 



By\t=0 — 



_( 1 -t- cos(107rz) z e [-0.1,0.1], 



elsewhere 



(A.; 



following iHood et al.1 l)2002a|) we expand By{C,t = 
0) in terms of Hermite polynomials By{^',t = 0) = 
J2'^=o'^rie~^ ^'^'^'TiniC /^) retaining only the fundamen- 
tal mode n = 0, i.e. to a fairly good accuracy initial con- 
ditions fEa. (|A.8p ) can be approximated By{^\t = 0) = 
aoe"^ /^'^ (note that TCq = 1). Therefore, substitut- 
ing the latter expression for By{^',t = 0) into the inte- 
gral (|A.4|I . and evaluating the integral analytically we ob- 
tain the solution of Eq. (|A.ip to the leading (fundamental 
71 = 0) order 



By = 



ao 



{z-CA{x)tf 



--[-TOt»]' ^ 

which coincides with the one obtained by iHood et al. I 
((2002a') . Here, cr = 0.033 and ao = l/b^a. For large 
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time, t, Ea. (|A.9p implies that the amplitude of the Alfvcn 
wave (e.g. amplitude of the By) decays as 

By = l[2nr^C'A^r/3]-'^'r'^^ (A.IO) 

i.e., in the case of a localized in z-coordinate Alfven pulse, 
Heyvaerts & Priest's e xponential decay is r eplaced by the 
power law By cx i"^/^ l|Hood et al. Il2n02a|^ . 



Although, the both derivations in the end lead to the 
same answer, in our opinion, the one presented here is 
more rigorous, concise and it gives the general solution 
fEa. ljA.4(l ^ and its asymptotic form fEa. ljA.5|l ) at once. 
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